Dynamical crystal creation with polar molecules or Rydberg atoms in optical lattices 
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We investigate the dynamical formation of crystalline states with systems of polar molecules or 
Rydberg atoms loaded into a deep optical lattice. External fields in these systems can be used to 
couple the atoms or molecules between two internal states: one that is weakly interacting and one 
that exhibits a strong dipole-dipole interaction. By appropriate time variation of the external fields, 
we show that it is possible to produce crystalline states of the strongly interacting states with high 
filling fractions chosen via the parameters of the coupling. We study the coherent dynamics of this 
process in one dimension (ID) using a modified form of the time-evolving block decimation (TEBD) 
algorithm, and obtain crystalline states for system sizes and parameters corresponding to realistic 
experimental configurations. For polar molecules these crystalline states will be long-lived, assisting 
in a characterization of the state via the measurement of correlation functions. We also show that as 
the coupling strength increases in the model, the crystalline order is broken. This is characterized in 
ID by a change in density-density correlation functions, which decay to a constant in the crystalline 
regime, but show different regions of exponential and algebraic decay for larger coupling strengths. 



I. INTRODUCTION AND OVERVIEW 

Recent experimental progress in controlled excitation 
of Rydberg atoms from cold gases [1-6] and in the pro- 
duction of ultracold polar molecules [7] has opened new 
opportunities for the study of systems with dipole-dipole 
interactions. Key features of these systems include the 
ability to control these interactions by varying the ex- 
ternal fields, as well as the strong dependence of the 
interaction strength on the internal state of the atoms 
or molecule. An important example of this dependence 
is the Rydberg blockade mechanism [8, 9] in which an 
atom in an excited Rydberg state can prevent nearby 
atoms from being excited to Rydberg states, because the 
strong interaction between the two excited atoms makes 
subsequent excitations non-resonant. This mechanism 
has been recently demonstrated in experiments [10, 11], 
and has important potential applications to the produc- 
tion of fast quantum gates with two or several neutral 
atoms [8, 9, 12, 13]. Moreover, the blockade produces 
rich and intriguing collective dynamics when atoms in 
a dense gas are excited to Rydberg states [1-6, 14-16]. 
This cooperative behaviour is manifest in a suppression 
of the proportion of excited atoms with increasing density 
(due to the presence of more atoms within the blockade 
radius), and also in an enhanced laser coupling to collec- 
tive Rydberg excitations. 

In the context of Rydberg atoms, there has been a lot 
of interest in the possible production of crystalline states, 
where a small fraction of atoms in excited Rydberg states 
arrange themselves in a regular array [17-20]. In partic- 
ular, when ground and excited states are coupled by a 
laser field, the resulting effective models correspond to 
crystalline states in the limit of weak coupling, with the 
fraction of excited atoms determined by a competition 
between the interaction energy and the detuning of the 
laser field from the atomic transition. 



Motivated by this experimental and theoretical pro- 
gress, we consider here a collection of either polar 
molecules or neutral atoms loaded into a deep optical 
lattice, e.g., with one particle per site. We assume that 
tunnelling in the lattice is negligible on the timescale of 
the experiment. Using dressed rotational states for the 
polar molecules or Rydberg states for the atoms we show 
how a microscopic model can be realized where two in- 
ternal states of the particles are accessed, one of which 
exhibits a large dipole-dipole interaction, whilst the other 
is weakly interacting. We show that this system allows 
the realisation of crystalline formations in the internal 
states, with a periodicity related to the lattice period by 
a rational fraction. The regularity of the system makes 
it possible to realize crystals with high filling factors of 
strongly interacting states. In addition, the use of po- 
lar molecules makes it possible for the crystalline states 
to be long-lived, aiding in the measurement of charac- 
teristic correlation functions. We show that in a ID 
tube, increasing the coupling strength between the inter- 
nal states at a fixed detuning will break the crystalline 
order. This is characterized in ID by a change in the 
density- density correlation (DDC) functions, which de- 
cay to a constant in the crystalline regime, but show 
different regions of exponential and algebraic decay for 
larger coupling strengths. 

We study in detail the dynamical preparation of states 
with crystalline order of different periodicities in the 
internal states by an appropriate time-variation of mi- 
crowave or laser couplings. As was shown in previous 
studies for Rydberg atoms in an optical lattice modelled 
by nearest-neighbour interactions [14-16], it is not possi- 
ble to produce crystalline states simply by switching on 
the coupling field. Instead, this is made possible using 
an adiabatic variation of the parameters in the effective 
model. Using an extended version of the Time-Evolving 
Block Decimation (TEBD) algorithm [21, 22], we calcu- 
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late the coherent time-dependent many-body dynamics 
of the excitations. We show that for system sizes and 
parameters corresponding to reahstic experimental con- 
figurations with either polar molecules or Rydberg atoms 
we obtain fidelities > 99% compared with the ideal tar- 
get states, i.e., the ground states of the effective model. 
Whilst the fidelity could be reduced by incoherent pro- 
cesses, e.g. decay of excited states, we estimate that these 
rates are small, and that the basic structure of the crys- 
talline state will be robust. Note that the time-dependent 
preparation of crystalline phases here follows the same 
philosophy as that for production of Rydberg crystals 
in reference [20]. In focussing on the case of atoms or 
molecules initially trapped in a lattice, we concentrate 
here on the case with a high proportion of particles in 
the strongly interacting case. This is in contrast to refer- 
ence [20] , which treats the case of a dilute gas of Rydberg 
excitations that is more relevant for a frozen Rydberg 
gas not trapped in a lattice. Our study of the correlation 
functions in the crystalline and non-crystalline phases is, 
however, also relevant for the case of [20]. 

The remainder of the paper is organized as follows: 
in section II, we introduce the effective long-range spin- 
model that describes both the Rydberg atoms and polar 
molecules in appropriate parameter regimes. Then we 
briefiy discuss in section III extensions to the TEBD al- 
gorithm that we use to study the model. In section IV, 
we analyse the ground states of this effective model and 
identify a target state of our time-dependent excitation 
process. The ground states are characterized as belong- 
ing to different quantum phases via the characteristic 
behaviour of DDC functions. In section V, we discuss 
a scheme to prepare crystalline states with high fidelity 
by an appropriate time-dependence of the external cou- 
pling field parameters. In section VI, details of two dif- 
ferent physical realizations of the effective model, with 
polar molecules (section VI A) and Rydberg atoms (sec- 
tion VI B) are provided. Finally, section VII, provides 
the conclusion and outlook. 



II. THE SYSTEM 

We consider particles loaded into an optical lattice in 
a regime where the lattice is sufficiently deep that tun- 
nelling can be neglected on the timescale of the exper- 
iment. We assume that the particles have two internal 
states, one of which exhibits a strong dipole-dipole inter- 
action, whilst the other is weakly interacting. This can 
be realized for realistic experimental parameters either 
with 

1. Polar molecules in an external electric field, where 
a microwave coupling of rotational states can give 
rize to dressed states with different dipole-dipole 
interactions, or 

2. Neutral atoms in an external electric field coupled 
via a laser to excited Rydberg states, which in con- 



trast to the ground states exhibit a strong dipole- 
dipole interaction. 

Below we summarize the key points of each of these im- 
plementations. Further details for polar molecules and 
neutral atoms coupled to excited Rydberg states are 
given in section VI A and section VI B respectively. 



A. Polar molecules 

Polar molecules in their electronic-vibrational ground- 
state manifold possess a number of features that make 
them appealing to design and implement specific spin- 
models. In addition to their electric structure, allowing 
e.g. for optical trapping, they have a rich internal (rota- 
tional) symmetry, with level spacings in the GHz regime, 
and negligible decoherence. They possess a permanent 
dipole-moment do (typically on the order of a few De- 
bye [23, 24]) along the internuclear axis, which makes it 
possible to strongly couple their rotational excitations via 
static and microwave fields and gives rise to dipole-dipole 
interactions between different molecules. 



FIG. 1. Schematic illustration of a system of polar molecules 
loaded into subsequent sites of an ID optical lattice with site 
separation a. The molecules are aligned by a static electric 
field in e^^-direction, Edc, and exposed to two time-dependent 
microwave fields (t) , E2 (t) , both linearly polarized along 
the z-axis. 

We consider a ID optical lattice (lattice spacing a) 
along the e^^ axis with one polar molecule trapped at 
each lattice site in the lowest vibrational state (j)k (x) (see 
figure 1), which can be realized either by loading an in- 
sulator state of polar molecules, or beginning deep in 
a Mott Insulator state with two atoms per site and as- 
sociating them to form polar molecules [25]. Without 
external fields, the rotational states of the molecules do 
not possess a net dipole moment. However, finite dipole 
moments along the z-axis in the rotational states can be 
induced by applying a static external electric field 
in that direction. The key idea of our implementation 
is that a suitable choice of E2{t) makes it possible to 
cancel the induced dipole moment in one of the dressed 
states, and we denote the complete state of an atom in 
this dressed state at site k as \a)k- The single molecule 
ground state, which obtains a strong dipole moment from 
the static electric field, is denoted as In the case 

that £^dc and E2{t) are properly tuned, two molecules at 
distinct sites k and / that are in the states and 
lead to an energy shift due to the dipole-dipole interac- 
tion potential of Vdd/|^ ~ ^P- section VIA we give ex- 
ample experimental parameters showing that an energy 
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shift Vdd ^ 0.7(io/47reoa^ ~ lOfikHz can be achieved. 
Note that we can therefore assume that the molecules 
remain in their motional ground state, since the band 
separation (trapping frequency) ujt is typically of the or- 
der of many tens of kHz up to 100 kHz in a realistic 
scenario. Finally, the states \P)k and \a)k on all sites k 
are coupled by a weak microwave field El{t) with Rabi 
frequency and detuning ^, which gives rise to the effec- 
tive microscopic many-body model that is described in 
section H C below. 



B. Neutral atoms coupled to Rydberg states 

Neutral atoms in an external electric field exhibit 
strong dipole-dipole interactions when in certain Ryd- 
berg states with relatively large principal quantum num- 
ber. For alkali atoms in states with principle quantum 
number n = 17, it is possible to obtain a dipole moment 
of more than one kiloDebye (see section VI B for details) . 
By finding isolated Rydberg states, the internal dynam- 
ics of each atom can be reduced to a two level system, 
where weakly interacting atoms in the ground state are 
coupled via a laser to these states with large dipole-dipole 
interactions. 

We assume that either atoms are loaded into a deep 
optical lattice with tight radial confinement, so that we 
have one atom per site along a ID tube, or that we have 
weak radial confinement and many atoms per site (as oc- 
curs, e.g., when an optical lattice is applied along only 
one direction). In each case, we can assume that the 
ground-state atoms are confined to a single modefunc- 
tion at each given lattice site, (/)/c(x). In the case of a 
single atom per site, this is the Wannier function for the 
lowest Bloch band, which is analogous to the case of the 
polar molecules described above. For many atoms per 
site, interactions can populate higher bands, and 0/c(x) 
will depend on the interactions. However, we assume that 
the timescale for excitation to the Rydberg state is much 
shorter than timescales corresponding to the atomic mo- 
tion, and that therefore we can neglect population of any 
modefunction other than the initial modefunction 0/e(x) 
for either the ground or Rydberg states. We then assume 
that we can have at most one excitation per lattice site, 
which in the case of large single-site populations is made 
possible by the Rydberg blockade mechanism, when we 
ensure that the Blockade radius is larger than the occu- 
pied region at any single lattice site (see section VI B). 
This makes it possible to describe the system with an 
effective spin- 1/2 model, assigning a state on each lattice 
site to the existence or non-existence of an atom in an 
excited Rydberg state on that site, labelled for site k by 
the states \I3)]^ and \a)k^ respectively. 

In comparison with polar molecules, the dipole-dipole 
interaction of the strongly interacting state for Rydberg 
atoms is much larger. If we choose a principal quantum 
number n = 14, then Vdd ^ TfiGHz (assuming a lattice 
spacing of a = 400 nm, see section VI B for more details). 



This corresponds to much faster experimental timescales. 
However, the polar molecules are much longer lived than 
atoms in excited Rydberg states, and open up different 
opportunities for experimental observation of the crys- 
talline states we discuss here. 



C. The effective model 

In either the polar molecule or neutral atom case we 
can derive the same effective microscopic model for the 
system in a frame rotating with the frequency of the 
field coupling the two internal states. We assume that at 
most one particle per site can be in the strongly interact- 
ing state, and represent the two possible states (weakly 
and strongly interacting) at each lattice site as two spin 
states. We label these states for site k as |ce)/c}, 
where \l3)k is the state exhibiting strong dipole-dipole 
interactions, and \a)k the non-interacting state. The 
Hamiltonian can then be formulated as (/i = 1) 



k k^m 
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with hk = \p){p\k being the local "number operator" 
for the interacting state at site k. The Pauli matri- 
ces are defined via = \P) {a\k -\- \ol) {P\k and = 
1/3) (/3 1 /c — \a){a\k^ where the latter is related to the num- 
ber operators via hk = {cf^ + 1)/2. Vt denotes the ef- 
fective Rabi frequency and 5 the detuning of the laser, 
which excites the particles, l^d is the dipole-dipole in- 
teraction energy shift between neighbouring sites. Note 
that we assume the particles to remain in the motional 
states ^/e(x), and also neglect corrections to the |A: — m|~^ 
decay due to the finite width of (j)]^ (x) (see section VI for 
more details). 

Below we will first study the ground state properties 
of this model, which will give us target states for the 
dynamical preparation process. We will then study the 
time dependence of this system starting from a many- 
body state where all of the particles are weakly interact- 
ing ((^^ \oL)k)i and investigate the adiabatic preparation 
of crystalline states as the parameters are varied. 

It is interesting to note that for an infinite system, the 
Hamiltonian can be written as 

Hs = nY,^l^\ (C(3)ydd -5)Y,^l 



Vc 



dd 



k 



^ \k- 



k^m 



my 



(2) 



where we have removed constant terms, and C(3) = 
Y.k^~^ ^ 1-202 denotes Apery's constant [27]. Note 
that in the case of ^ = C(3)Vdd, the system is equivalent 
to an Ising model in a transverse field with dipolar in- 
teractions. Ground states of this model have been also 
studied in [281. 



4 



III. NUMERICAL METHOD 

To calculate the many-body dynamics and ground 
states of Hamiltonian (1), we make use of the TEBD al- 
gorithm as introduced in [21, 22]. This algorithm makes 
possible the near-exact integration of the Schrodinger 
equation for ID lattice and spin Hamiltonians based on 
a matrix product state ansatz, and has been applied to 
a range of such Hamiltonians with next-neighbour in- 
teractions. There are also applications of matrix product 
state algorithms to dissipative systems [29-31], and these 
methods have been incorporated within density matrix 
renormalization group algorithms [32, 33]. There is also 
a strong effort to extend these ideas to higher dimensions 
[34]. 

The modification in our work is the extension of these 
methods to the treatments of finite-range interactions. 
For the purposes of this study, we include long-range in- 
teractions up to an arbitrary distance of / sites, by imple- 
menting a routine for swapping site indices in the matrix 
product state representation. This adds a computational 
cost of 0{l) basic operations compared to the standard 
TEBD calculation. Here we work with finite systems of 
size A/", and always properly represent interactions over 
the full length of the system, i.e. / = A^ — 1. This method 
of swapping also allows us to consider periodic boundary 
conditions for the dipole-dipole interactions. Note that 
we always perform convergence tests in the matrix prod- 
uct state bond-dimension x and other numerical param- 
eters [35]. For larger system sizes, swapping sites can be- 
come inefficient, in part due to the higher required values 
of X- To implement interaction terms over substantially 
larger distances than in this study, more efficient algo- 
rithms should be used, e.g., algorithms involving use of 
matrix product operators [36-38]. 



IV. GROUND STATES OF THE EFFECTIVE 
MODEL 

We begin by studying the key ground state properties 
of the effective model of Hamiltonian (1). This will allow 
to us to determine in which parameter regimes crystalline 
order will appear, and to characterize the states that will 
be the target states of the dynamical preparation process 
discussed next in section V. 

The Hamiltonian we study here has strong similari- 
ties to another model treated in several recent studies, in 
which the dipole-dipole interactions were replaced with 
van der Waals interactions arising for certain Rydberg 
systems [16-20]. In [17] it was shown that for Q = 
that system reduces to a classical spin model exhibiting 
a second order phase transition from a paramagnetic to 
a crystalline phase at (5 = 0, and we observe similar be- 
haviour here. Due to the dipolar long range character 
we expect crystalline states with different periodicities of 
the excitations to the strongly interacting state to appear 
for (5 > and small values of Q. The detuning plays the 
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FIG. 2. A sketch of the key features of the phase diagram for 
the system with Hamiltonian (1). For ^ < and ^/Vdd ^ 
1, the ground state involves all of the atoms in the weakly 
interacting state in order to minimize the interaction energy 
(marked by the blue lobe). As the detuning ^ > is increased, 
it becomes favourable to add excitations to the system, but 
this competes with an increasing interaction energy. This 
competition gives rise to crystalline phases with periodicities 
determined by this competition (marked by the brown lobes). 
For large the crystalline order is broken, but some 

particles exist in excited states (white area). 



role of a chemical potential for strongly interacting states 
(or external magnetic field in the spin model), and it be- 
comes favourable to excite states with a dipole-dipole 
interaction on the lattice. However, the long-range in- 
teractions will compete energetically with the detuning, 
leading to crystal periodicities that decrease with increas- 
ing 5/Vdd- Below we show that as the coupling fl/Vdd is 
increased, the crystalline order is broken, in favour of a 
paramagnetic phase in the spin language, with no long 
range density-density ordering. These key features of the 
phase diagram are sketched in figure 2. 

We compute the ground state of equation (1) using the 
TEBD algorithm in several different parameter regimes. 
As a first step, we look at the total number of particles 
in a strongly interacting state 



N 



(3) 



In figure 3 we present a plot of Nf^ for several detunings 
S and effective Rabi frequencies We plot N^s on a grid 
with spacings AS = 0.1251/dd and AQ = 0.0251/dd. 

As expected, for small 0./Vdd there is essentially no 
occupation of the strongly interacting states in regions 
with negative detuning, because a state with no \f3)k ex- 
citations is energetically favoured. In contrast, a positive 
detuning leads to a reduction of the total energy when 
occupations of strongly interacting states are added to 
the system, and in regions with large positive detun- 
ing we find states with occupations of \P)k on all lat- 
tice sites. In the region between we observe varying 
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FIG. 3. Occupation number of strongly interacting states, 
Nj3^ for ground states of the effective model in a 50-site sys- 
tem. The Hamiltonian is given through equation (1) and open 
boundary conditions are assumed. Results are shown for sev- 
eral values of detuning 5 and Rabi frequency and plateaus 
of zero, half, and full filling of the lattice are clearly visible for 
^ < O.lVkd (see text). The contours of = 23 and = 27 
are drawn as lines on the plane. 



excitation numbers depending on the choice of param- 
eters. In particular, there is a large plateau in the region 
0-5 Vdd ^ ^ ^ 2]/dd and Vt < O.lVdd^ in which Nj^ remains 
close to N/2 = 25. In the plot the contours of fillings 
slightly smaller {Nj^ = 23) and larger {Nf^ = 27) than 
N/2 are drawn as lines on the 5-Q plane. Note that the 
N/2 plateau lobe is centred around a value slightly larger 
than (5 = 1, which is approximately the point where the 
second term in equation (2) vanishes and our system be- 
comes an effective Ising spin model. 

The half filling plateau disappears with increasing Q 
and we find a linear interpolation between regions with 
zero and full filling in this regime. This is consistent 
with a breakdown of the crystalline structure, and we 
note that in the limit 1] ^ oo for finite 5 we ex- 
pect the filling to be constant at N/2 since the state 

(8)/c i\f^)k ~ \(^)k) /V2 becomes the ground state of the 
system. 

Note that the grid is too coarse-grained to resolve 
plateaus of smaller fractions of the filling. However, we 
observe several kinks in the particle number at small val- 
ues of We show below, e.g., that the feature at 
S ~ 0.25Vdd corresponds to a crystalline state with filling 
N/3. This will verify the lobe-like structure as depicted 
in figure 2, however we find that the plateau for filling 
A^/3 is significantly smaller, at values of 1] lO^^'Kid- 
For values of S ^ 2.25Vdd we also observe some kinks, 
which correspond to regular patterns in the lattice with 
two subsequent sites being in the strongly interacting 
state, followed by one site being in the state \a)k- 

In order to illustrate the claim that the plateaux seem 
to belong to crystalline states, we can investigate the 
fluctuation in the total occupation number of the strongly 
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FIG. 4. Fluctuation of the total occupation number of the 
strongly interacting state in a 50 site system with open bound- 
ary conditions as a function of for 5 — l.lVdd- The inset 
shows a zoom into the region 0.05Vdd < ^ < 0.5Vdd. 

interacting state, 

N 

AN^ = Y,{m)-Nl (4) 

k,l 

In figure 4 we show for S = 1.1 Vdd and 0.05Vdd < 

Q < 5Vdd- We find that AA/"^ decreases significantly with 
decreasing Q and a crossover behaviour is centred around 
ft ^ 0.2 Vdd- For Q < 0.2 di crystalline state with nearly 
no fluctuations AN 13 < 3 <C = 50 is present. For large 
ANi^ converges to a value close to N/A = 12.5, which 
is the correct result for finite 5 and 1^ ^ 00, i.e. where 
®k ~ \(^)k) /a/2 becomes the exact ground state 

of the system. Note that the quantity is related 

to the Mandel Q-factor Q = AN^/N^ — 1, which has 
been measured in recent experiments for the number of 
Rydberg excitations in small ensembles of cold atoms [39] 
and quantifies the deviation from a Poissonian number 
statistic [40]. 

To study ground states in more detail for different pa- 
rameter regimes, and demonstrate clearly the crystalline 
behaviour we evaluate DDC functions, defined as 

^['^ = {hihij^k) - {rii){ni+k)' (5) 

In a finite system with open boundary conditions, the 
exact values of Q^^ will differ from the site i where the 
DDC is evaluated. However, the decay behaviour of the 
correlation functions for large k should be independent 
of the site index z in a sufficiently large system. To dis- 
tinguish ground state phases by their decay behaviour, 
we will first calculate ground states of a system with pe- 
riodic boundary conditions, in which the DDC becomes 
site independent. We then drop the site index and write 

Qk- 

In figure 5(a) we show a characteristic DDC in the 
parameter regime where the filling of the lattice is ap- 
proximately A^/2, choosing a detuning of ^ = l.ll^d. 
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FIG. 5. DDC functions Qk for the ground states of Hamilto- 
nian (1) for 6 = l.lFdd and (a) Q = O.lKid, 0.2ydd, and (b) 
Q = O.SVdd, 0.4Vdd- The insets show the modulus \Gk\ on a 
logarithmic scale. Here we show results for a system with 50 
sites and periodic boundary conditions. We find states with a 
marked crystalline structure of periodicity 2 in panel (a) and a 
phase transition at 0.2Vdd < < O.SVdd in terms of the decay 
behaviour of In contrast to panel (a), in (b) \Gk\ decays 
exponentially. In the case of Q = 0.4Vdd the long-range decay 
behaviour {k > 10) follows a power-law decay oc Ck~^ (with 
C ^1.7x 10"^). 



with both n = O.lVdd and n = 0.2ydd. The inset in fig- 
ure 5(a) shows the modulus of the correlation function, 
\Qk \ on a logarithmic scale. In the regime of the half filling 
plateau we indeed find crystalline structure with a pro- 
nounced zigzag pattern. The long-range crystalline order 
is indicated by the decay of \Qk\ to a non-zero value as a 
function of the site index k. We find that for increasing 
N the value of the constant to which \Qk \ decays becomes 
independent of the system size. 

We note, already from the comparison in figure 5(a) 
that as ^/Vdd increases the strength of the long range 
crystalline order decreases. If we continue to increase 
^/Vdd5 we find that this order is broken entirely, as shown 
in figure 5(b), with the same parameters, but for Rabi 
frequencies Q = O.Sl^d and Q = 0.4Vdd- In these cases 
we also observe peaks with periodicity 2, however, these 



occur only for very small values of /c < 5. As k increases, 
\Qk \ decays to zero, indicating the breakdown of the long- 
range order. This appears as a phase transition from 
crystalline order to essentially a "paramagnetic phase", 
in which the correlations typically decay exponentially. 
Here we note a special feature of the correlation func- 
tions due to the long-range interactions, however, which 
is most clearly visible in the results for Q = OAVdd in 
figure 5(b). There we observe that the long-range or- 
der has been replaced by exponential decay at short dis- 
tances {k < 10), but that this changes from exponential 
to a power law-decay. By fitting a function we find that 
in this region \Qk\ ~ C/k^^ with a constant C. This is 
the same power law decay as the interaction term, as 
is consistent with similar results derived for other long- 
range spin models [28, 41]. Physically, the correlations 
initially decay exponentially due to the primary influ- 
ence on a given particle being from its nearest neighbours 
(analogous to exponentially decaying correlations in the 
paramagnetic phase of an Ising model). However, once 
correlations arising from nearest neighbour interactions 
become sufficiently small, long-distance correlations will 
be dominated by effects of the direct interaction between 
distant particles. Thus, long range correlations are ex- 
pected to decay at large distances with the same power 
as the interaction [41]. 

We find similar behaviour for crystalline states at other 
densities, including occupation periodicities larger than 
two. When the detuning is decreased from 6 ^ O.SVdd 
it becomes less favourable to occupy \l3)k states in the 
system and therefore in the case of small crystalline 
phases with larger spacing between these occupations ap- 
pear in order to reduce the dipole-dipole interaction en- 
ergy. For example, in figure 6 we show the DDC indicat- 
ing an occupation of the states \P)k on every third site 
for Q ~ 10~^Vdd and S = 0.25Vdd- To make this period- 
icity clearer, we present results for a system with N = 60 
sites. As in figure 5, in figure 6 we observe the charac- 
teristic long range order in the DDC, combined with a 
marked peak on every third lattice site. Again we find 
a phase transition away from crystalline order, with an 
exponential decay of \Qk\ exhibited for larger values of 
ft. Note that the crystalline phase with 1/3 filling occurs 
in a much smaller region of parameter space than the 
crystalline phase with 1/2 filling. 

In what follows, we are interested in time-dependent 
creation of these crystalline states in the case of polar 
molecules or Rydberg atoms loaded into an optical lat- 
tice. These systems will typically occupy ~ 50 lat- 
tice sites, with open boundary conditions generated by 
the unoccupied ends of the chain. Thus, rather than at- 
tempting to derive the full phase diagram for an infinite 
system, we focus on the case of open boundary conditions 
relevant to the experiments. Although the correlation 
functions are modified by the boundary conditions, we 
show below that the characteristic signatures of different 
phases can still be clearly identified. 

In figure 7 we show the DDC evaluated from site num- 
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FIG. 6. DDC functions Qk for the ground states of Hamilto- 
nian (1) for 6 = 0.25ydd and (a) Q = O.OlFdd, 0.02ydd, and 
(b) Q = 0.04ydd, 0.05Vkd. The insets show the modulus |&| 
on a logarithmic scale. Here we show results for a system 
with 60 sites and periodic boundary conditions. We find a 
marked crystalline structure in panel (a) with periodicity 3 
and a phase transition at around Q ~ O.OSVdd indicated by 
the exponential decay behaviour of \Gk\ in panel (b). 
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FIG. 7. DDC functions evaluated from site index i 
50-site system with open boundary conditions, Q]^K 
suits shown are for the same parameter regime as in figure 5. 
Despite strong boundary effects, a crystalline phase is still 
distinguishable from a phase with such structure via the be- 
haviour over short ranges /c < 10. In the crystalline phase the 
long-range decay follows a power law behaviour proportional 
to Ck~^ {C ^ 1.7 X 10~^ is a constant). 



axis, and shows that the boundary depicted corresponds 
to a phase transition. 

Note that the DDC functions can be measured ex- 
perimentally both for Rydberg atoms and for polar 
molecules. Detection of these correlation functions is dis- 
cussed in more detail in section VI. 



V. DYNAMICAL CREATION OF 
CRYSTALLINE STATES 

A. Time-dependent state preparation 



ber 5 of a A/" = 50 site system, i.e. Q^. ^ with open 
boundary conditions. We choose the same parameters 
as were shown in figures 5(a) and (b), including crys- 
talline phases with periodicity 2. Despite obvious effects 
of the open boundary conditions we can clearly distin- 
guish crystalline behaviour = O.lFdd, O-^'Kid) from the 
case were exponentially decaying density-density corre- 
lations are present {Q = O.S'Kid, 0.4l^d)- While in the 
crystalline phase we find very little decay for increas- 
ing k until the effects of the boundary become apparent 
{k > 20), in the case of > 0.3 the correlations rapidly 
decay below values of 10~^. As in figure 5 we still observe 
a long-range power law decay for k > 10 proportional to 
k-\ 

In relation to figure 2, the results presented in figures 5- 
7 give evidence for the crystalline phases at filling fac- 
tors 1/2 and 1/3 on the lattice, and show for each that 
there exist transitions as Q is increased to non-crystalline 
phases. This justifies that the lobes are finite along the Q 



We now discuss the dynamical creation of crystalline 
states in our system. The natural initial state for the ex- 
periments corresponds to having all of the particles in the 
non-interacting states \a)k^ which would be the ground 
states for neutral atoms, or a single dressed rotational 
state for polar molecules. In terms of the effective model, 
this state is the ground state of the Hamiltonian (1) for 
large negative detuning and small Q. The question is now 
whether one can time-dependently vary the laser param- 
eters ft and (5, i.e. find a trajectory in the 1^-5-plane, 
so that this initial state is adiabatically transferred to a 
crystalline state of the form discussed in section IV. 

In order to guide our discussion of this question, we 
first consider a small model system consisting of 8 sites, 
where exact diagonalization of the Hamiltonian (1) is 
straightforward. We can then estimate possible paths for 
the adiabatic transfer based on the size of the energy gap 
AE between the ground state and excited states. Dur- 
ing the adiabatic ramp, we would like to remain as 
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FIG. 8. Shaded plot of the energy gap /\E between the ground 
and the first excited state as a function of and 5. The 
results are from exact calculations for an 8-site system with 
Hamiltonian (1). The lines indicate the contours of constant 
A path, which consists of three segments and yields an 
adiabatic passage to a crystalline ground state, is sketched as 
orange arrow (see text). 



large as possible to suppress non-adiabatic transitions to 
excited states of the effective model. In figure 8 we show 
a shaded plot of the energy gap between the ground state 
and the first excited state as a function of Vt and J, with 
lines of constant energy gap marked in the plot. Close 
to = we find a large region where the energy gap 
assumes small values and where one can recognize lobe- 
like structures. Especially visible are lobes at values of 
5 ^ 2.25Vdd- As mentioned before these correspond to 
states with regular patterns of two subsquent sites being 
in the strongly interacting state. In contrast to the lobe 
for filling fraction N/3 of the strongly interacting state at 
S ^ 0.25Fdd7 these are enhanced due to the larger increase 
in interaction energy when adding single occupations of 
the strongly interacting state to the system. Therefore, 
these lobes are more robust when increasing Q. In gen- 
eral, the system is too small to observe lobes of filling 
fractions and smaller and these are only visible as 
small kinks for S ~ 0.25Vdd- With increasing Rabi fre- 
quency ft the energy gap increases. The requirement of 
having at any instance of time a sufficiently large energy 
gap is most easily fulfilled if one chooses to "steer" a path 
in parameter space around those regions where the gap 
is the smallest. A simple path that seems to achieve this 
is, for example, to first increase Q at constant large neg- 
ative detuning S and then increase at a constant large 
Q. Both of these two steps can be achieved by changing 
the parameters S and Q rapidly, due to the large gap 
present at all times. Afterwards, ft can be decreased to 
a small value on a line of constant S. This path consist- 
ing of these three segments is marked on figure 8 by an 
orange arrow. 

Now we would like to test these insights by investigat- 
ing the time-dependent propagation of the system along 
this path in parameter space. In order to investigate the 
scaling with the system size, we study systems of size 
20, 30, 40 and 50 lattice sites. We note that the large 



sizes here are already typical of likely experimental se- 
tups. As an example, we choose to prepare a crystalline 
state with excitation periodicity 2, which is the ground 
state of the our model for the parameters Q = O.OSVdd 
and 5 = l.lVdd- We perform real time simulations for a 
path consisting of three segments, as sketched in figure 8. 
We start ai Q = O.OSVdd and S = —3Vdd- There, the in- 
ner product between the ground state of the Hamiltonian 
and the state in which all atoms are in the \a)k state is 
sufficiently large for all system sizes. The segments of the 
path are: 

1. We increase ft from 0.05Vdd to 0.5Vdd at J = — 3 
with a variation rate Aft/ At — 0.05F^^^. 

2. We increase ^ from — 3Vdd to l.fKid at ^7 = O.S'Kid 
with a variation rate A5 / At = 0.05V^^. 

3. We decrease ft from O.S'Kid to O.OS'Kid with several 
variation rates AQ/At. 

Due to the large gaps in the energy spectrum it is pos- 
sible for segments 1 and 2 to obtain states with a high 
fidelity F = |(^g(^, 1 in the comparison 

of our evolved state, with the ground state of 

the effective model with the corresponding parameters, 
1^^(1^,(5)) on that path. The most sensitive part of this 
process is the third segment. In figure 9 we show the 
fidelity F of the adiabatically evolved state along seg- 
ment 3 as a function of ft{t) and for various variation 
rates AQ/At. On = 30 lattice sites, we find that for 
the three rates Aft/ At = 5 x 10"^l/dd, 1 x lO-^Kid, and 
5x 10~^Vdd5 this fidelity drops below values of 80%, which 
occurs when crossing the phase transition from figure 5, 
i.e. in the regime 0.2Vdd < ft < O.SVdd- At this point 
the ground state gap becomes too small and excitations 
to higher excited states analogous to Landau-Zener tun- 
nelling processes occur. However when using a rate of 
Aft/ At = 1 X lO~^Vdd we find that the fraction of the 
state that is lost into higher excitations remains reason- 
ably small and the final fidelity with the crystalline state 
at ft = O.OSVdd is larger than F ^ 99%. Thus, we have 
shown that a crystalline state can be produced with very 
high fidelity in a finite system of 30 sites. The timescale 
which is required for all three segments is of the order 
of 500/Fdd7 and we will show that this is experimentally 
realisable for our two physical implementations below. 

In an experiment with either polar molecules or neu- 
tral atoms coupled to excited Rydberg states in an optical 
lattice, the final state could be characterized by measur- 
ing DDC functions. These can be most easily detected 
directly by in situ imaging experiments [42] . In the case 
of polar molecules, the crystalline structure will be vis- 
ible in noise-correlation measurements of state-selective 
momentum distributions of molecules released from the 
lattice [43-45]. In these correlation measurements, the 
periodicity of the crystalline structure would be present 
as a peak at the corresponding reciprocal lattice vector. 
In the case of Rydberg excitations, direct measurement 
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FIG. 9. The fidelity of tiie adiabatically evolved state with 
the current ground state, when changing Q from 0.50Vdd to 
0.05Vdd at constant S = 1.1 Vdd in segment 3 of the adiabatic 
preparation path. Results are shown for different variation 
rates AQ/At within a 30-site system with Hamiltonian (1). 
For decreasing variation rates the final fidelity for the crys- 
talline state (at Q = 0.05 Vdd) increases towards 1. 



FIG. 10. The fidelity of the time-evolved state with the target 
crystalline state, at = 0.05 Vdd, and S — 1.1 Vdd as a function 
of the time required for the segment 3. Here we show results 
for several system sizes of 20, 30, 40 and 50 lattice sites. This 
fidelity decreases with the increasing system size and with 
decreasing preparation time. 



of DDC functions could also be made by imaging the 
sample on a channel plate detector [20]. 

An important question is: to what extent the same 
fidelities we have observed for a system size of 30 lat- 
tice sites can also be achieved for larger system sizes. In 
general we expect the energy gap between ground and 
excited states to become smaller with increasing system 
size, and we expect the measure of the fidelity to become 
more sensitive due to the exponentially increasing size of 
the Hilbert space. To analyse how much this affects the 
final fidelity, in figure 10 we plot the fidelity achieved af- 
ter ramping the laser parameters through all segments of 
our preparation procedure, as a function of the time re- 
quired for the final segment 3. We consider system sizes 
of 20,30,40 and 50 sites. 

We find that for all system sizes the fidelity is > 94% 
when choosing a variation rate slower than AVt/At = 
1 X 10~^F^^^, which corresponds to an operation time of 
approximately 450 /Vdd on segment 3. In the case of an 
even slower rate of AVt/ At = 5 x 10~^V^^^ which requires 
a time of approximately 900/'Kid7 the fidelity becomes 
larger than 98%, even in the case of a system with 50 
lattice sites. This fidelity is very high, especially given 
that the Hilbert space dimension is 2^^ in the latter case. 



B. Comparison with experimental parameters 

Let us compare the time required to reach a final crys- 
talline phase with experimentally accessible timescales, 
which are in general limited by the finite lifetime of 
the strongly interacting state. In the scheme presented 
above, at the slowest variation rate of AVt/At = 5 x 
10~^l^^^ on segment 3, a total operation time of the or- 



der lOOO/'Kid is required. Note that this time can be 
further dramatically decreased by choosing shorter paths 
or by not utilising a constant ramping rate but rather re- 
ducing it with time while approaching the critical point 
of the phase transition. 

In the case of the polar molecule implementation (see 
section VIA for more details), typical parameters give 
Vdd ^ lOfikHz, assuming a permanent dipole moment 
(io = 1 Debye and a lattice spacing of a = 400 nm. This 
translates to a time for the complete preparation of ap- 
proximately 0.1s, which is significantly shorter than re- 
cently measured lifetimes of molecules in an optical lat- 
tice of 8 s [46]. By choosing polar molecules with larger 
dipole moments like LiCs (do ^ 5.5 Debye [24]), the 
preparation time can be further reduced down to the or- 
der of a few ms. 

In the case of the Rydberg atom implementation (see 
section VI B for more details), if we choose a principal 
quantum number n = 14, then Vdd ^ ThGHz (assum- 
ing a lattice spacing of a = 400 nm). Hence, the to- 
tal operation time corresponds to approximately 140 ns 
in the experiment. This value is well below the typical 
lifetime of a n = 14 Rydberg level of Lithium, which 
can be found to be approximately 2.2 /is [47]. Note that 
these estimates neglect the effects of blackbody radia- 
tion, which can redistribute population amongst excited 
levels. For Rydberg atoms, in general the lifetime in- 
creases with increasing n and furthermore the required 
operation time decreases as n~^. However, the disad- 
vantage of large principal quantum numbers is the small 
level distances Ainter and Aintra (see figure 13), which 
both decrease with n. This limits the experimentally ac- 
cessible values of 6. In the case n = 14 (a = 400 nm) 
at an electric field strength of F^i = lOOkV/m, we can 
estimate by using the results known from the hydrogen 
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atom Ainter ^ 4.6 X lO^V^d and Aintra ~ QOKid. The 
field strength is thereby well below the Inglis- Teller limit, 
which is in this case 319kV/m [48]. This justifies the va- 
lidity of our two- level model, derived in section VI B, in 
this regime. 

This comparison to the lifetime of a single Rydberg 
atom is oversimplified, as in the presence of ex- 
cited atoms, the rate of single decay events will increase 
roughly proportional to Nf^ [49]. However, as the system 
size (and hence A^^) becomes larger, the effect of a sin- 
gle decay on the crystalline structure should be reduced. 
Whilst the many-body state fidelity will be very sensi- 
tive to a single decay, DDC functions should be robust, 
at least on length scales smaller than those separating 
places where spontaneous emissions have occurred. The 
effect on the many-body state of single decay events also 
depends on when during the excitation process they oc- 
cur. A full treatment of this excitation process including 
decay events could be achieved by studying the dynam- 
ics of a master equation including spontaneous emission 
events. 

Therefore and from the results in figure 10 we conclude 
that the adiabatic passage should be experimentally real- 
isable for both our implementations of the effective spin- 
model with a size of ~ 50 sites. For larger system sizes, 
longer timescales will be required to reach as high fideli- 
ties as we have obtained here. However, the properties 
of the final states will typically be determined in an ex- 
periment by measurements of the correlation functions. 
In contrast to the state fidelity, correlation functions (to- 
gether with the physical character of the state) become 
more robust to small non-adiabaticities (especially those 
resulting in localized defects) in larger systems. 



VI. PHYSICAL IMPLEMENTATIONS OF A 
LONG-RANGE SPIN MODEL 

A. Polar molecules in an optical lattice 

In the following we describe how to implement Hamil- 
tonian (1) with a system of polar molecules. We con- 
sider N molecules confined to the lowest motional level 
on subsequent sites of a deep ID optical lattice, with 
modefunction ^^(x) on site k. We assume that the lat- 
tice is sufficiently deep that tunnelling between wells of 
the lattice is strongly suppressed on typical experimen- 
tal timescales. In particular we consider molecules in 
their electronic vibrational ground-state manifold with a 
closed shell structure ^S, which possess permanent dipole 
moments (e.g. "^^K^^Rb like in recent experiments [7]). 
The setup is depicted in figure 1. We use a strong and a 
weak microwave field linearly polarized along the z-axis 
and an additional static field aligned along the same di- 
rection. Given the fact that we consider samples of the 
order of a few to hundred /im one can neglect the spa- 
tial variation of the respective static and microwave field 
along the optical lattice axis, which is aligned in the e^^- 



direction. 

The low energy spectrum of a single molecule at site k 
in our setup is well described by rotational excitations, 
governed by a rigid rotor Hamiltonian coupled to the 
fields via an electric dipole interaction [50-53] 

HUt) = BJl - dl (Ede + El{t) + El{t)) . (6) 

Here, B is the rotational constant typically of the or- 
der of tens of GHz {h = 1) [26], is the angular 
momentum operator, the z-component of the dipole 
operator and E^^^ ^i(^) = exp(— icjit) + c.c, and 
E2{t) = E2 ex.p{—i(jj2t) +C.C. the static, weak and strong 
external microwave fields, respectively. 

The eigenstates of are denoted | J, mj)k with eigen- 
values J{J -\- 1) {h = 1), and mj = —J,..., J being 
the eigenvalues of the projection onto the z-axis, J^. 
The dipole operator dl couples rotational states with dif- 
ferences A J = ±1, and we recall that the eigenstates 
of equation (6) have no net dipole moment. The non- 
vanishing elements of the operator d^. are given by ( J zb 
l,mj\ dl\J,mj)k = do{J,mj; 1 , 0| J=b 1, mj)( J, 0; 1,0|J =b 
1,0)V'(2J + 1)/(2(J±1) + 1) in terms of Clebsch-Gor- 
dan coefficients (Ji, mi; J2, J, m), where do denotes 
the permanent electric dipole moment along the internu- 
clear axis, which is typically of the order of a few Debyes 
[23, 24]. 

By applying electric fields one lifts the degeneracy of 
the state manifolds with J > and induces a finite dipole 
moment in the dressed rotational states. Since the fields 
are aligned along the quantisation axis, mj is conserved, 
and we focus on states with mj = in our setup, which 
are connected to the ground state. The induced dipole 
moments of two molecules at distinct sites k and / then 
lead to an energy shift due to the dipole-dipole interac- 
tion potential. 

In order to describe this, we need to consider many 
molecules on the lattice. Molecules in an internal state 
denoted a can then be represented by field operators 
^^cr(x), obeying the standard bosonic or fermionic com- 
mutation relations. We consider a situation where no 
motional modes other than 0^ (x) for each site k are pop- 
ulated during the experiment, which is valid provided 
that the band separation ut is substantially larger than 
the laser coupling parameters and interaction energies. 
Thus, we can expand the field operators as 

^.(x) = ^0,(x)^^ (7) 
k 

and identify individual molecules occupying each mode 
described by the operators g^. Since all rotational states 
we consider here have mj = 0, they are only coupled 
by an interaction proportional to d^df, which can be 
integrated over the modefunctions in each lattice site. 
The resulting many-body Hamiltonian of a system for 
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FIG. 11. Level scheme for a system of polar molecules in 
a combination of static and microwave fields, cf. figure 1. 
(a) Energy levels Ej of a molecule in a static electric field 
Edc- The field orients the dipoles, i.e. induces finite per- 
manent dipole moments in each rotational state |J)dc- (b) 
A strong linearly-polarized microwave field E2 (t) couples the 
states |l)dc and |2)dc the Rabi- frequency Q2 and detuning 
A2. Fine-tuning the microwave- field (with respect to the 
static field) makes it possible to cancel the permanent dipole- 
moment in one of the two dressed states, and cf. 

da ^ 0. A second (weaker) microwave field Ei(t) finally cou- 
ples the ground state \/3) = |0)dc to the dressed state with 
vanishing dipole \a) with an effective Rabi- frequency Qi and 
detuning Ai. 



making the rotating wave approximation leads to the two 
dressed states \a±)k = cos(^/2)|l)dc — sin(6>/2)|2)dc with 
the dressing angle = arctan(— 21^2/^2) (0 < ^ < 27r), 
where the state \a-)k (|a+)fe) corresponds to ^ < tt 
{6 > tt). The corresponding energies of the dressed states 
are = Ei - A2/2 ± v'A^ +41^2/2. 

Finally, we couple the state \a)k = \(^-)k and the 
state = |0)dc by the weak electric field with Rabi 
frequency fti = {a\dl\P)kEi and detuning Ai = ui — 
{Ea — Eq). Thereby we assume that the field El{t) is 
much weaker than £^2(^)5 ^^^t the states \a) and 
are unaffected from Ef{t). Going to the rotating frame 
and making the rotating wave approximation the many- 
body Hamiltonian becomes time-independent and reads 
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(9) 



witho-j; = |/?)(a|fc + |a)(/?|fc and Ufe = |/?)(/?|fc. In the last 
term of equation (9), Vd is the dipole-dipole interaction 
matrix in the rotating frame, which in the two-molecule 
basis {\a)k\a)i, \a)k\/3)i, \(3)k\a)i, \/3)k\l3)i} is 
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molecules can be written as 
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(1) V- 
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with lattice spacing a, and where V^^^^ ~ l/(87reoa^). 
Note that we neglected corrections to the shape of the in- 
teractions, which are of the order 0(/h/(<^|^ — ^|)), where 
/h denotes the characteristic harmonic oscillator length 
for the optical potential well on each site. 

A suitable choice of external electric fields then makes 
it possible to strongly couple/dress locally two (effective) 
rotational levels, such that only one level has a strong in- 
duced dipole moment, while at the same time the transi- 
tion dipole moment remains negligibly small. This results 
in dipole-dipole interactions in equation (8) reducing to 
the diagonal form of equation. (1). 

We now detail our appropriate choice of the electric 
fields to obtain the desired dressed states. The key idea 
is to apply the static field E^c to induce a finite dipole 
moment in the rotational states, leading to the dressed 
states |J)dc, defined by [BJl - dlE^c]\J)dc = Ej\J)^c 
(see figure 11(b), and to then couple the two states |l)dc 
and |2)dc by a strong microwave field £^2(^)5 ^^^h that the 
dipole moment vanishes in one of the resulting dressed 
states. The field £^1 (^) tuned to the transition between 
the two states with Rabi frequency 1^2 = (2|<^||l)dc£^2 
and detuning A2 = 0J2 — (^2 — ^1) (see figure 11). Di- 
agonalizing the corresponding two-level Hamiltonian and 



Here, the non- vanishing elements of Vd can be evaluated 
in the rotating wave approximation from the bare single- 
particle dipole-operator elements dij = (i|<i||j)dc- 



= {cos\0/2)dii^sm\0/2)d22) 

+ 2sin^(l9/2)cos^(l9/2)d?2 
= doo (cos2(^/2)dn + sin' (^72)^22) 
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Our goal is to find optimal values of the microwave 
dressing angle and static field strength £^dc so that 
the desired Hamiltonian [equation (1)] is realized with as 
small an error as possible. Though from equation (13) we 
see that there is no parameter choice where is exactly 
zero, we can find regimes where the model is implemented 
up to very small errors. At this point we note that sub- 
tracting a multiple of the identity from Vd amounts to 
an overall shift of the zero of energy, and thus we can 
equivalently write 
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FIG. 12. Matrix elements of the dipole-dipole interaction ma- 
trix Vd in static and microwave fields, Edc and E2{t). (a) 
Interaction strength for the strongly interacting state, 6v (in 
units of do), (b) The ratio \v^a ~ '^a^l/^'^ diagonal 
element of Vd and (c) of the off-diagonal element, v'^^/dv, 
as a function of the static- field strength do Edc /B and the 
microwave- dressing angle ^ on a logarithmic scale, (d) The 
relative deviation, 6Vd/Sv, of the induced dipole-dipole cou- 
pling with respect to the model interaction as a function of 
do Edc /B and ^ on a logarithmic scale. The solid lines show 
contours of constant value. 



In Figure 12 we plot the matrix elements Sv^ \v^^ — 
f^^l, and v^^ as a function of the microwave field angle 
and the electric field strengths up to E^c = 20do/B. In 
figure 12(a) we show the magnitude of the interaction en- 
ergy for the strongly interacting state, Sv. Figure 12(b), 
and figure 12(c) show the ratio of the other interaction 
matrix elements to Sv, i.e. \v^^ ~^a^\/^^^ l^a^l/^^- 
In panel (d) we plot the (normalized) total deviation from 
an exact implementation of the Hamiltonian [equation 
(1)]. This is quantified as SVd/Sv, given by 

5v^ = \\v;-5v\p}kmimk{m 

= \/te-<^)^+2(t;f;)2. (16) 

We find that for field strengths ~ ^^B/do, fine tun- 
ing the microwave field at 6 ^ O.dir makes it possible 
to obtain deviations SVd ~ 10~^Sv, while the interaction 
strength for the strongly interacting state is Sv ~ 0.65^0 
(see contour line in panel (a) ). Thus within small errors 
we obtain the model interaction with a typical timescale 

of Vdd = '^^^'^dd ^ Sv/4:7reoa^ ~ lOfikHz for typical pa- 
rameters (a = 400 nm, do = O.SDebye). An alternative 
approach using circularly polarized coupling fields is also 
possible, which can give rise to an implementation of the 
model for appropriately tuned parameters [54]. 



B. Neutral atoms coupled to Rydberg states 

In this section we show how to implement spin Hamil- 
tonian (1) by making use of neutral atoms coupled to 
excited Rydberg states. We consider alkali atoms con- 
fined in an optical lattice in a single tube, with inter-site 
distance a oriented along the x-direction. Perpendicular 
to the lattice a homogeneous electric field, F = F^ie^ is 
applied, together with a laser which couples the ground 
state 1^) of each atom to a specific chosen Rydberg-Stark 
state. 

i. Dipole-Dipole Interaction 

While the atomic ground state is only slightly shifted 
by the external electric field, the level structure of high- 
lying Rydberg states is strongly altered compared to the 
field-free case. Nevertheless, the azimuthal symmetry of 
the atomic Hamiltonian is preserved such that m (the 
quantum number of the z-component of the electronic 
orbital angular momentum L^) remains a good quantum 
number [48]. States with large quantum defect, i.e. usu- 
ally s- and p-states, sustain a second-order Stark shift 
while the degenerate states with higher angular momen- 
tum exhibit a linear Stark effect, reminiscent of the hy- 
drogen atom [55]. The derivative of the energy with re- 
spect to the electric field strength directly yields the elec- 
tric dipole moment of the given state. While the precise 
level structure and therefore the dipole moment of the 
Stark states depends on the element under consideration 
the exact results known from the hydrogen atom usually 
provide good estimates for the quantities of interest. 




-► 



FIG. 13. Sketch of the level structure of an alkali metal atom 
as a function of the electric field strength. The energy of 
hydrogen-like states grows linearly with the field strength Fei. 
Low angular momentum states with quantum defect exhibit a 
second order Stark effect. In our scheme we couple the ground 
state 1^) of an atom to the Rydberg state \dn), which exhibits 
the strongest Stark shift within the given n-manifold. The 
magnetic field strength and the laser parameters are set such 

that Aintra < Ainter, Aintra ^ \S\ and Aintra > ^0- 

A sketch of a Rydberg-Stark level structure is shown in 
figure 13. For sufficiently small field strengths, i.e. below 
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the Inghs-Teher hmit where the atom becomes ionized, 
the levels are grouped in manifolds which can be labelled 
by the principal quantum number n. In our study we 
focus on the energetically highest state of such a manifold 
whose azimuthal quantum number is m = 0. Employing 
the formula for hydrogen, the electric dipole moment of 
this state, which we denote as \dn), evaluates to 



dn = dne^ = -eaon{n - l)ez, 



(17) 



with e and ao being the charge of the electron and the 
Bohr radius, respectively. 

Two atoms excited to the state \dn)j. and \dn)i at dis- 
tinct sites of the lattice k and /, are expected to interact 
strongly via the dipole-dipole interaction potential, as al- 
ready for n = 17 equation (17) gives a dipole moment of 
more than one kiloDebye. We intend to work in a regime 
in which the interaction can be treated within first order 
perturbation carried out in single particle product states. 
To this end it is important to note that the dipole-dipole 
interaction commutes with + , and therefore only 
couples atomic pair states for which the azimuthal quan- 
tum numbers obey mi + m2 = m[ + m2 [8] . We assume 
that all of these couplings are far off-resonant, which can 
be assured by properly tuning Aintra and Ainter using 
the electric field. Moreover, accidental resonances can 
be avoided by restricting oneself to low principal quan- 
tum numbers n since here the two-particle energy spec- 
trum possesses large gaps which make near-resonant ex- 
citations unlikely. Employing -L e^c the interaction 
Hamiltonian between two atoms being in the product 
state \dn)]^ \dn)i can then be written as 



rjint 
^kl 



2V^ 



(1) 



dl 



dd 



\k-i\' 



\dn)k\dn)l{dn\k{dn\ 



(18) 



with V^^^'^ ^ l/(87reoa^), when the confinement length in 
each lattice site is much smaller than the lattice spacing. 

We now turn to the description of the interaction be- 
tween the laser and the atom. We employ a two-level 
approximation assuming that the ground state \g)j^ is 
coupled only to the state \dn)j.- This is justified as long 
as the modulus of the laser detuning S and the Rabi fre- 
quency r^o are much smaller than the energy gaps Aintra 
and Ainter- Making the rotating wave approximation, the 
interaction of the laser with an atom located on the k-th 
site can be written as {h = 1) 



^0 {\dn)k{9\k + \9)k{d„\k) - S\dn)k{d 



"n\k 



(19) 



with fto and S being the Rabi frequency and the detuning, 
respectively. 



2. Many- 



model 



We now consider the many-body physics when many 
atoms in an optical lattice are simultaneously coupled to 



the state |(in)/c- We assume that either atoms are loaded 
into a optical lattice with tight radial confinement, so 
that we have one atom per site along a ID tube, or that 
we have weak radial confinement and many atoms per 
site. In each case, we can assume that the ground-state 
atoms are confined to a single modefunction at each given 
lattice site, (/)/c(x). In the case of a single atom per site, 
this is the Wannier function for the lowest Bloch band, 
and in the case of many (weakly-interacting) atoms it 
is the solution to the Gross-Pitaevskii equation in the 
mean- field limit. 

In order to formulate the Hamiltonian of the many- 
body system, we introduce field operators ^^(x) and 

(x) , corresponding respectively to atoms in the ground 
and Rydberg states, and obeying the standard bosonic 
(or fermionic) commutation relations. We assume that 
the timescale for excitation to the Rydberg state is much 
shorter than timescales corresponding to the atomic mo- 
tion, and that therefore we can neglect population of any 
modefunction other than the (x) for either the ground 
or Rydberg states. In this case, we can approximately 
expand the field operators in terms of mode operators 
corresponding to each of these modes. 



k 

^^(x) = ^0fe(x)ffc, 



(20) 
(21) 



where both the initial field operators and the mode oper- 
ators, gl and f^, obey the appropriate bosonic commuta- 
tion or fermionic anti-commutation relations, depending 
on which atoms are used in the experiment. We can then 
write the Hamiltonian for a lattice with TV sites as 

N N 

Hn ^^0^ ^gkrl + givk - S'^hk 
k k 

N N 



2 f-f \k-l\^ 



(22) 



with fik = T^Vk^ Vdd ~ (i^/47reoa^. In the last term 
we have again neglected corrections to the shape of the 
interaction potential, which are of the order 0{li{/ {a\k — 
/|)), with the characteristic harmonic oscillator length 
for the deep optical potential well on each site. 

The quantity V^^ is relevant only in the case that we 
initially have many atoms per site, and gives the energy 
offset due to the dipole-dipole interaction if two Rydberg 
atoms are excited on the same site. We will assume that 
V^^ > Vdd/IA: ioik and l^o, {N5) < V^^, so 
that the occupation of a single site by two Rydberg atoms 
is far off-resonant and therefore the excitation of such 
a state is highly improbable. This phenomenon, which 
is known as the Rydberg blockade, makes it possible to 
restrict the maximum number of Rydberg atoms per site 
to one, such that hy^ can only assume the two eigenvalues 
and 1. 
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Secondly, to arrive at a Hamiltonian governing solely 
the dynamics of the Rydberg excitations, we eliminate 
the creation and annihilation operators gj^ and gk in equa- 
tion (22). In the case where we have an initial ground 
state with a single atom per site, this is simply a rela- 
belling of the state with a single ground state atom per 
site as the vacuum state of the operators Vk- In the case 
where the mode function at each site is occupied by the 
same macroscopic number A^g 1, on the other hand, we 
make the replacement gk^gl -> in which case we 

define the effective Rabi frequency as = Qoy^^N^. The 
model we derive is exact under the stated assumptions 
when we have a fixed number of atoms per site (e.g., be- 
ginning deep in the Mott Insulator regime for bosons or 
in a band insulator for Fermions), which is the ideal ex- 
perimental realization. An approximation to the model 
can also be obtained in a weakly interacting gas, but with 
some inhomogeneity in Q for different lattice sites, both 
due to position-dependence of the mean density (due to 
a trapping potential), and also due to number fluctua- 
tions between sites. In production of crystalline states 
via the adiabatic ramping process discussed below, we 
expect that the ramp is somewhat robust against these 
fluctuations, especially for high filling factor crystalline 
states. However, this depends on the details of the real- 
ization and the inhomogeneity in Q. 

Thus, we arrive at the spin Hamiltonian describing ex- 
citation of Rydberg states, 

N N AT . . 

^r.a = 0E-^^ + ^E-^ + T^E^' (23) 

k k k^^l 

where the Pauli matrix = + rk has been defined. 
This Hamiltonian is equivalent to our model system from 
equation (1), with the states corresponding to \a)k = \g)k 
and \P)k = \dn)k' 

The scheme presented in section V achieves a crys- 
talline phase of Rydberg excitations on very short time- 
scales. However, whilst the preparation can be very fast 
compared to the lifetime of the Rydberg states, this life- 
time will pose an ultimate limit to the stability of the 
many-body state that is produced. One way around 
this is to transfer the excited Rydberg states to alter- 
native stable ground states. However, an important is- 
sue arises in this case because of the momentum kicks 
that arise due to the interactions between excited atoms, 
which can create very significant excitations in the mo- 
tional states of the atoms. This is especially an issue 
towards the end of the chain, where the atoms will ex- 
perience a force primarily in the direction away from the 
centre of the chain. These kicks can be estimated clas- 
sically from the force acting due to the dipole-dipole in- 
teraction potential between excited Rydberg state. For 
crystalline states one can estimate that the net momen- 
tum kick during the crystal preparation time decreases 
as kji (X l/(/z)^, where i denotes the i-th excited Ryd- 
berg atom from the boundary and / the periodicity of the 



crystal, i.e. with N/f filling. Quantitatively, transitions 
to higher bound states are given by the size of matrix 
elements J Wn{x)wo{x)e'^^^^dx^ where {n > 0) de- 
notes the higher band Wannier mo defunct ions. Thus the 
transition effects become small in the case that kRa^, 
where ao denotes the oscillator length, is small. The 
oscillator length is typically of the order of tens of nm 
and we can estimate the momentum kick of the bound- 
ary excitation for a N/2 filling crystal as roughly 1/nm. 
Thus, this will be a problem for atoms at the end of the 
chain in the N/2 case, but decreases rapidly for atoms 
away from the edge, or a larger periodicity in the crys- 
tal. One way to circumvent this problem for the case 
where all lattice sites are singly occupied would be to use 
induced dipoles rather than direct coupling to Rydberg 
levels from the outset, as discussed in reference [56]. The 
idea is to make use of two stable ground states, one of 
which is coupled off-resonantly to a Rydberg level, pro- 
ducing an induced dipole with a strength depending on 
the parameters of the coupling and the static dipole of 
the Rydberg level. In this way, the coupling between 
ground and Rydberg states in our model is replaced by 
the coupling between two stable ground states. This can 
be achieved, e.g., via a two-photon transition. This re- 
sults in much smaller values of Q and 6 (and hence also 
the interaction strength between excited atoms), which 
could even be made comparable to or less than the band 
gap in the optical lattice. The atoms would then always 
remain trapped in the lattice. The disadvantage of this 
method would be substantially increased times for the 
total adiabatic process, however, the ratio of the prepa- 
ration time to spontaneous emission lifetime of the state 
could always be made large. 



This method would have the added benefit of a much 
smaller linewidth of the coupling field for the transi- 
tion between interacting and non-interacting states. Es- 
pecially for crystal periodicities other than two lattice 
sites, the narrow regions of detuning S/V^d for which the 
phases are found implies a strong condition on the sta- 
bility of the laser detuning in an experiment (and hence 
on the linewidth of the exciting laser). For example, a 
crystalline state with an excitation on every fourth site 
can be found in the region of ^/Vdd = 4 x 10~^ for 
0.06 < 5/Vdd < 0.09. Using Vdd ~ IhGUz from the 
previous example, this means that the range of detuning 
values in which states of the character is found is approx- 
imately 200 MHz, and we require that the laser linewidth 
is smaller than this value. This is realistic, however for 
crystalline phases with even larger spacings between the 
excitations, the requirements become even stronger. If 
we use dressed dipoles, then we would be using a two- 
photon transition to couple between the states, for which 
the detuning can be controlled very accurately. 
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VII. CONCLUSION AND OUTLOOK 



In conclusion, using adiabatic sweeps of microwave or 
laser parameters, it is possible to form crystalline phases 
of strongly interacting states in either polar molecules 
or neutral atoms coupled to Rydberg states in an opti- 
cal lattice. The effective model Hamiltonian exhibits a 
rich ground state phase diagram, with the crystalline or- 
der being destroyed as the coupling between states is in- 
creased. For typical experimental size scales and param- 
eters we showed that preparation of crystalline structures 
can be achieved on reasonable timescales, and the struc- 
tures can be detected by measurement of characteristic 
density- density correlation functions. 

The results we have obtained here motivate further 
study of the phase diagram with its lobe-like structure of 
different crystalline phases (as indicated in figure 8). This 
also opens the opportunity to extend ID time-dependent 
simulations to the study of crystal formation in regimes 
of significantly larger spacing than is studied here (e.g.. 



using matrix product operator methods [29, 30, 36-38]) 
and also to probe time-dependent transitions between the 
different phases. This will also allow connections to cur- 
rent work in which Rydberg excitations of randomly dis- 
tributed atoms are studied [20]. 
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